On numerical integration of coupled Korteweg-de Vries System 
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Abstract 

We introduce a numerical method for general coupled Korteweg-de Vries systems. The scheme is 
valid for solving Cauchy problems for arbitrary number of equations with arbitrary constant coeffi- 
cients. The numerical scheme takes its legality by proving its stability and convergence which gives 
the conditions and the appropriate choice of the grid sizes. The method is applied to Hirota-Satsuma 
(HS) system and compared with its known explicit solution investigating the influence of initial con- 
ditions and grid sizes on accuracy. We also illustrate the method to show the effects of constants with 
a transition to the non-integrable case, 



§1 Introduction 

^Coupled Korteweg-de Vries system equations form a class of important nonlinear evolution systems. Its 
m ^importance comes (physically) from the wide application field it covers and (mathematically) from includ- 
^ing both (weak) nonlinearity and third order derivatives (weak dispersion). It describes the interactions 
^o!f long waves with different dispersion relations. Namely it is connected with most types of long waves 
with weak dispersion (co(k) — > 0, (k) — > 0), e.g. internal, acoustic and planetary waves in geophysical 
hydrodynamics. 

It was introduced by A.Maxworthy, L.Redekopp and P.Weldman M in studying the nonlinear atmo- 
sphere Rossby waves. R.Hirota and J.Satsuma || give single- and two-soliton solutions to some version of 
the system. R.Dodd and A.Fordy || found an L-A pair for Hirota-Satsuma equations. S.B. Leble derived 
the cKdV system for different hydrodynamical systems with explicit expressions for the nonlinear and dis- 
persion constants ||. He also developed the approach to the cKdV integration. S.B. Leble, S.P.Kshevetskii 
||, |5| used the system in investigation of nonlinear internal gravity waves. A. Perelemova |J used it in 
description of interaction of acoustic waves with opposite directions of propagation in liquids with bubbles. 

Others deal with integrability of the system (Dodd R. and Fordy A. H) from a Lax pair point of view, 
S.B. Leble in Walquist-Estabrook theory. Foursov MV || described a new method for constructing 
integrable system of differential equations that reduced to cKdV equations. Oevel W. |{J considers inte- 
grable system of cKdV and found an infinite hierarchy of commuting symmetries and conservation laws in 
involution. Zharkov A Yu. jnj obtained a new class of integrable KdV-like systems. Metin Gurses, Ata- 



lay Karasu [11 1 found infinitely many coupled system of KdV type equations which are integrable. They 



also give recursion operators. In studying the Painleve test classification of the system, Ayse (Kalkanli) 



Karasu [12| found new KdV systems that are completely integrable in the sense of WTC paper. He was 



looking for the integrable subclass of KdV systems given by Svinolupov [TSl]. The later has introduced 



a class of integrable multicomponent KdV equations associated with Jordan algebras. John Weiss [14 



derived the associated "modified" equations for HS system and from these the Lax pair is also derived. 



B.A.Kupershmidt |15[ showed that a dispersive system describing a vector multiplet interacting with the 



KdV field is a member of a bi-Hamiltonian integrable hierarchy. 

The significant achievement in numerical solution of the single KdV equation starts from the famous 



paper of Norman Zabusky and Martin Kruskal |TB[. It develops the idea of soliton solutions set for the 



integrable equations and enlighten the problems of effective integration scheme elaborating. The paper 
launched the numerous investigations and inventions in this field. Perhaps, the last publication that 
develop applications of recent theoretical achievements in numerical integration schemes is based on the 
notion of isospectral deformations |18| . Recently a multisymplectic twelve points scheme was produced 



17j . This scheme is equivalent to the multisymplectic Preissmann scheme and is applied to solitary waves 
over long time interval. 

Shaohong Zhu |]TU[ had produced a difference scheme for the periodic initial-boundary problem of the 
Coupled KdV (Hirota- Satsuma case) system. He use the inner product of the discrete function to obtain a 
scheme keeping two conserved quantities. His scheme is a nonlinear algebraic system for which a catch-ran 
iterative method is designed to solve it. 

The coupled KdV system representing most possible physical application ( related to the weak nonlinear 
dispersion ) to be considered in this work takes the following general form 

{0 n ) t + c„ {9 n ) x + ^ g m kn Ok (0 m ) x + d n {0 n ) xxx = , n,m, k = 1,2,3, N, (1) 

k,m 

where 9 n (x, t) is the amplitude of the wave mode as a function of space x and time t respectively. The 
constants c n are the linear velocities and g m kn, e n are the nonlinear and dispersion coefficients. 

In the present work we introduce a numerical tool for solving coupled KdV system which is a develop- 
ment of the two step three time levels as Lax- Wendroff scheme |5|, ^(J . Proving the theorem about stability 
and convergence of the scheme gives the opportunity to use it for different applications like Cauchy prob- 
lems for arbitrary number of equations and a wide class of initial conditions 8 n (x, 0). We consider in our 
problem an infinite domain while the initial condition goes quickly enough to zero following the relation 

1 + \x\)\6(x, 0)\dx = < oo 

keeping in mind the choice of smooth and integrable function. As an important corollary of the theorem 
one obtains conditions that have to be taken in account in choosing grid sizes. This numerical method is 



checked by applying it to HS system for which a good number of explicit solutions exist |2T[ . We examine 
also the effects of equations coefficients and conditions of the problem on the solution. 

In the section 2 we introduce the difference scheme for arbitrary number of coupled KdV equations. 
We investigate stability and prove the convergence giving the condition have to be taken in account in 
choosing the grid sizes and how they are related. In section 3 we analyze HS system with two-parameters 
one-soliton explicit solution. The numerical method is applied to HS system and compared with the 
explicit solution. We analyzed the effects of the two parameters and initial condition on the form of the 
resulting solitons as will as on accuracy and show the results by figures. We also produced (numerically) 
a multi-soliton solution for HS system and used the conservation law to estimate the expected number 
of solitons which agreed that we already obtained. Proving stability and convergence besides testing the 
results for HS system allows us in section 4 to use scheme for general cKdV system. Hence we illustrate 
by plots the results of applying the scheme to slightly nonintegrable cKdV systems and other for a system 
with non-smooth initial conditions. 



2 The numerical method 



2.1 The difference scheme 

For the cKdV system (1) we introduce a numerical ( finite-difference ) method of solution. A scheme 
which is two steps three time levels similar to the Lax-Wendroff one ||, ^0). The usual Lax-Wendroff is 
modified such that the order of the first derivative becomes of order 0(Ax A ). The approximation of the 
nonlinear terms is changed in such a manner that the integral of 9 2 be a conserved one. The approach 
gives a solution that can be considered as some generalized solution, in the sense of Shwartz distribution 
theory, where the dispersion constants vanishes. This scheme is suitable to a nonlinear equations and is 
valid for n equations with arbitrary coefficients. The scheme can be simply derived beginning from Taylor 
series expansion as 

(9 n )i +1 = (e n )i + At((e n ) t y i + o[(Atn (2) 

where i and j are used to locate a point in the discrete domain and At is the time step while the subscript 
t means time derivative. Substituting for [6 n ) t in (2) using the system (1) to obtain 

(0 n )l +1 = (6 n )i - At (c n (6 n ) x + J2 9mkJ k (0 m ) x + d n (e n )\ + O [(A t) 2 ] . (3) 

\ k,m / j 

The difference scheme is elaborated applying Lax idea for a half time step and leap frog method to the 
remaining half time step. In both steps {6 n ) x and {0 n ) xxx are replaced by forth order O (Ax 4 ) and second 
order accurate O (Ax 2 ) central difference expressions. Hence (3) gives the following difference scheme 

((e n )i H - (dn)i) l\ + c n ((9 n y i+1 - (e^Q /2h + J> mfc7l (e k )i ({e m y i+1 - (e^Q /2h 

k,m 

+ e„ (ie n y i+2 - 2 (9 n )i +1 + 2 (9^1, - (^ n )|_ 2 ) /2h 3 = 0,e n = (d n - c n h 2 /6) , (4.1) 

where n, m, k are the modes numbers; i and j are discrete space and time respectively. The time step At is 
replaced for simplicity by t while h denotes spatial step. The equation (4.1) is accompanied with discrete 
equation for the intermediate layer as: 



(>n)f 1 - (0„)J) It + c n ([6 n )\ + J - (9 n )li) /2h + Y,9mkn * {{oj^ - [dj^) /2h 

k,m 

+ e n {[OniXl ~ 2 infill + 2 (0„)£? - (e n )\ + _l) /2h 3 = (4.2) 
2.2 Stability and convergence analysis 

For simplicity of the analysis we start by considering one equation of the system and give the details of 
stability and convergence. Then we apply the idea to the general cKdV system because it is rather close 
to that for one KdV equation but more bulky. 

2.2.1 Stability analysis for KdV scheme 

Consider one KdV equation of the system (1) 



e t + c9 x + g ee x + de xxx = o 



(5) 



Note again that the investigation we perform can be generalized for the case of any finite number of modes. 
Considering the numerical scheme applied for the equation (5) 

(ot 1 - 01) It + C (ej +1 - eU)/2h + 9 m + , - oL)/2h + e « 2 - 2e{ +1 + 29U - e{_ 2 )/2h 3 = o (6) 

Let us select a suitable norm. For this multiply equation equation (5) by 9 and integrate to yield 

1 d ro ° 



2dt 



9 A dx = or 



f 

J — c 



9 dx = const, 



hence, by definition of L 2 norm,(||#|| 2 ) 2 = 9 2 dx ,ft may be written as 
norm ||0|L is conserved and the equation can be treated in the L 2 norm. 



,) 2 = const., i.e. the 



Now we will prove stability with respect to small perturbations ( because we consider nonlinear equa- 
tions ) of initial conditions. Strictly speaking it is the boundness of the discrete solution in terms of small 
perturbation of the initial data. So let us consider the differential 



d9j +1 = J2(d9l +1 /d9l)d9l , 



+ 1, ... 



(7) 



' d01_ 2 



for equation (5) denoting 7V+ 1 = {d9j +1 /d9 3 r } , d9{ 



d9\ 
d9\ 
d9l 
d9i 



> and use = £ {d9lf h 



+i 



+2 f 



Rewrite also the relation (7) in the matrix form 

d9{ +1 = Tl+ X d9{ = T^TW -1 = Y[T r d9°, 

r 

where d9° is a small perturbation of initial data and the subscripts are omitted for simplicity. 



Stability requires the boundedness of the product Yl T r in a sense that the norm 



is bounded 



by some constant, i.e. 



< C . Here C is a constant, and the matrix norm is a spectral norm. 



For this the sufficient condition is ||T r || < e aT where a is a constant, that is independent of r. The 
case ||T r || < e a ^-< h )* T is a sufficient condition of stability also, but only if |a(r, h)\ < const < oo. If 
|a(r, h)\ < const, including r, h — > 0, for some dependence r = f(h), then we can say about conditional 
stability. Namely this kind of stability will be climbed below. 

To calculate T r , rewrite the scheme (6) in the form 9j +1 = 9 j + 1 (9 j i+2 , 9 J i+1 , 9j, 9\_ x , 9\_ 2 ). 

So 



(Ti +1 ) ir = 5 hr - (cr/2h) [6 i+1 , r - 6^] - (gr/2h) \9{ (6 i+1>r - <W) + 6 i>r {9{ +l - 

- (er/2/i 3 ) [5 i+2 ,r - 25 i+ i jr + 2Si- 1>r - 5i^ 2 ,r 

Rewriting (8) in terms of the identity (E), symmetric (S) and anti-symmetric (A) matrices 
Ti +1 = C + S j+1 + Ai +1 

{S 3+1 } hr = -f£ {{el - 9l +1 ) 5 t+1 , - {91 - 9U) 8^ + 28 hr \9{ +1 - flf_J) 



(8) 



i + ou) Si.!,) 



er 



2h 3 



[$i+2,r — 2S i+ i : r + 25j_ lr . — <5j_ 2 ,i 



||^ +1 || < \g\ r max(|^| , |^|) ,0i, = - Of] /(h) and 9% = [9f +1 - 6ff_J /(2/i) 
one arrive at 

ii am ii - M T I/171 l c l r 3 lei r 

11 11 ~ h i 1 1 /i /i 3 



| T i+i|| 2 = || ptf+i)* = || - + [E + + 

< 1 + 2||^ +1 || + (||^ +1 || + ||^ +1 ||) 2 



II II Q I I \ 2 

\g\ |c| 3 |e| 



< 1 + 2 |g| r max |#^| + r 2 I 2 |g| max + ^ max | + ^ + ,, 

2 Igl max \9l A + — max I + — H — — - ) 

i 1 21,11 /i i 1 11 /i /i 3 / 



which is a necessary condition of stability. The scheme is stable if a < constant < 00 in spite of r, h — > 0. 
This is a conditional stability of the scheme. It means that it is required for stability that r — > more 
faster than /i — * 0, or 

r < (constant), h 6 , constant < 00 (9) 
Therefore, for small enough r we can simplify the expression for a 

x,i I 



2|#| max r +r 



3e 
X 3 



In practical calculations the time step r should be chosen so that it would satisfy r (f§) 2 *^o = 0(1), where 
to is the time of simulation (0 < t < t ). In future, when we shall be suggesting some better numerical 
scheme, we will essentially use our observation that stability depends only on the dispersion terms. And 
now we will try to accomplish our short investigation of the scheme (6) by strong proving of the numerical 
scheme convergence. 



2.2.2 The proof of the KdV scheme convergence 

Now we prove that a solution of equation (6) converges to a solution of (5), if the exact solution is a 
continuously-differentiable one. Let us denote by 9(x,t) a solution of the equation (5). We substitute 9\ = 
9(xi,tj) + v{ into (6), v\ is a error between the difference solution 9\ and the exact solution 9(x i) tj). We 
obtain the equation for v\ 

{vl +1 - vj) It + c (vj +1 - vU) /2h + g9(x u tj) (v{ +l - «f_J /2h + gvj (9(x l+l) t j ) - 9(x l _ 1 ,t J )) /2h 

+ gvl - ^-1) /2h + e (vj +2 - 2vj +1 + 2vU - v{_ 2 ) /2h 3 = 
- ((9( Xl , t j+1 ) - 9(x t , tj)) It + c (9(x l+u t j ) - 9(x i _ 1 ,t j )) /2h + g9(x l} tj) (9(x i+lj t j ) - 9(x l „ 1 ,t j )) /2h 

+e (9(x l+2 , tj) - 29(x l+1 ,t j ) + 29(x i _ l ,t j ) - 9(x t _ 2 , tj)) /2h 3 ) 

Let us take into account that 

vl - r {c(vi +1 - vl 1 )/2h + g9(x h tj)(v{ +1 - vl 1 )/2h + gvi(9(x i+1 ,tj) - 9(x^, t J ))/2h 

+e(vi +2 - 2vj +1 + 2vU - vl 2 )/2h*) = £ {^%v> k 

k 



Using the operator T^ +1 introduced above, this equation may be rewritten in the form 



vt 1 - 



E ( Tj+1 )iA h + gvl (vUi ~ V U l 2h = 



k 



- ({0(xi, t j+1 ) - 9(xit j+1 )) /r + c (6(x i+1 ,tj) - eixi-^tj)) /2h + gO(xi, tj) (e(x i+1 ,tj) - Oixi-^tj)) /2h 

+e (6(x i+2 , tj) - 28(x i+1 ,t j ) + 28(x i _ 1 ,t j ) - 6(x t _ 2 , tj)) /2h 3 ) . 

The right part of this relation is a quantity of order 0{r + h 2 ). So, we can write 

W +1 " E ( T ' + % 4) h + 9vl (v{ +1 - vU) /2h = 0(r + h\ 

v k J 



or 



i+i 



E - rf!, fl = 9vl Vi+1 2 ^ U - 0(r + tf). (10) 

k 



We finally arrive at the inequality that compare the norms: 

ill ^ \9\ IL..7II 2 , , u2 



< ^\\v J \\ +0(r + h 2 ) 
To explain how this estimate was obtained, follow the expressions 



2 



x 2 \ 2 

I"* -1 I** 



2\ 



11/1 = (E (//) 2 *J * l«l IE (^^^TT^J "I + 0(T + " 

< Isl (E N) 2 * * E W) 2 M ' 3 + 0(r + ft 2 ). (11) 

Using the Schwartz inequality ||AB|| < \\A\\ \\B\\, the formulas (10) could be transformed as 

||^+i|| < ||T i+1 || H^'ll +r \\f j || < ||T j+1 || ||T J '|| ll^- 1 !! +r(||T i+1 || \\f j + ||f ||) < 

|| T i+i|| || T i|| || T i-i|| p'-2|| +r (|| T j+i|| || T j|| ||ji-2|| + || T i+i|| + n^'ii < 

e arJ |K|| +r(e ar(j - 1) ||/°|| + e ar( ^ 2) H/'H + ••• ||/ J ||) < 
e ari \\v°\\ + Mmax(\\f k \\) < 

e^||t; ||+M^|4||^ +1 |r + O(r + / i 2 )^ , M = \2 ~\ - ( 12 ) 

To derive (12), we have used the iteration of the first of formula (we substituted the formula into itself, 

but for index less than 1) and using (11). Then we have utilized the formula for a sum of geometric series. 
Farther the inequality obtained in (12) has a solution 

< ^1 - ^1 - 4m|4 (e ar \\v°\\ + MO(t + h 2 )) j / ( 2M j4) 



If we take into account (9), and use||u°|| = 0,we obtain 
\\v j+1 \\ < I 1- 



4M M »^/ \ , f 2M \g\ 



IMO(t + W) / - 
/is I \ h 2 



( 1 _^ 0(T + fc . ) )) / (2^2l) = M . 0(T + fc . 



The constant M is bounded, in spite of j — > oo, because of jr < oo. Therefore, the convergence is proved. 
2.2.3 The coupled KdV scheme 

The numerical scheme for the system of the cKdV (4.1)-(4.2) is also conditionally stable and convergent 
one. The proof for this scheme is close to that given before, but a bit bulky. We deal with a vector: 

u={e 1 e 2 ... e N Y 

as a dependent variable instead of the simple variable 9 in the case of one KdV. For this vector case the 
norm used has the form 

/ n \ § 

\ i=i 

The conditions connecting time and space steps for this scheme look also similar but with different constants 

81 * r 3 

max(|e n |) * - — — *t = 0(1). 

n 4 * Ai iz 

3 Checking the numerical method 

The numerical method is tested by applying it to Hirota-Satsuma system. Namely the two parameters 
one-soliton explicit solution is used [^1[ . 

3.1 Analytic solution ( explicit formula ) of Hirota-Satsuma system 

Darboux transformation (DT) that account a deep reduction for this specific HS case of cKdV is used 



21] to obtain explicit solutions to HS system. The Lax representation of the HS equations is based on 



the matrix 2x2 spectral problem of the second order. For this problem the deep reduction scheme [[Hj] is 



applied (with the help of the conserved bilinear - forms) and supports the constrains on the potential while 
the iterated DT are performed. The iterated DT in determinant form and the covariance of the bilinear 
forms with respect to DT under restrictions gives N soliton solution of HS system. The system of HS we 
use to check the scheme has the form: 

(6 1 ) t - 0.25 (6 1 ) 3x - 1.5 (6 1 ) x fa) + 3 (6 2 ) x (6 2 ) = (13.a) 

(e a ) t +0.5(e 2 ) ta + 1.5(e a )„(e 1 ) =0. (13.b) 
This system has two-parameters one soliton solution: 

0i = -2m 2 (-l + d 2 + 2dSin{\ 1 ) * Sinh{X 2 ))/{dCos{\ 1 ) + Cosh(X 2 )) 2 , 

9 2 = (2 + 2d 2 y 5 m 2 /(dCos(X 1 ) + Cosh(X 2 )), (14) 
Ai = .5m 3 t + mx and X 2 = .5m 3 t — mx. 



with real constants m, d. For small \d\ this solution is a smooth function but for \d\ > 1 poles appear .The 
following figures show some choice of m and d to show the effect of these two parameters on the solution. 
Figure 1 show that, for constant d, the amplitude is proportional to m while the wave width is inversely 
proportional to it. Figure 2 show that, for the given m, d affects on soliton shape, namely the first mode, 
while the amplitude of the second is inversely proportional to d. 




The first mode (u) 

m=5, d=0 

- - - m=3, d=0 



The second mode (v) 

m=5, d=0 

- - - m=3, d=0 




Fig. 1. For a constant d, the amplitude is proportional to m while the wave width is inversely 

proportional to m. 



The second mode (v) 
m= 5, d= 0.5 
m= 5, d=0 




Fig. 2. For the same m, d affects on soliton shape, namely the first mode, while the amplitude 

of the second mode is inversely proportional to d. 



3.2 Calculations by numerical scheme and comparison results 

HS system (13) is solved numerically using the scheme (4) with initial condition from formula (14) (t = 0) 
and the results are compared with the explicit solution. It is found that, keeping the restriction on the 
choice of r and h and relation between them, the initial wave modes amplitude affects on the accuracy 
of the results . Also the error decreases as the mesh is refined. Namely smaller amplitude ( of order one ) 
gives better results as shown below in figure (3). It gives the percentage error calculated as follow 



'lError 



\Explicit solution — Numerical solution\ 



Initial amplitude 

We relate the error to the initial amplitude to show the physical significance of the error. 



First mode, %Error, t=2 sec 
Ao = 2 , Ao=3.4 



= 5 — 




4 




3 j 




2 - 












, 1 





x (x=0 is the peak point) 



Second mode, " < > Error. t=2 sec 
Ao=2, Ao=3.4 



-5 O 5 

x (x-0 is the peak point) 



Fig. 3 % error is proportional to the amplitude(A). 



The plots show that the error is proportional to the amplitude (A), where as shown in figure the maximum 
relative error in the case (A=2) is 1 % while in the case (A=3.4) is 3 %. The reason may be due to the 



higher velocity in the larger one, hence more interactions impact. It also shows that the error increases near 
the peak points. The reason of these oscillations in plots appearance is that the numerical and analytical 
plots intersect over the space domain. 



4 Applying the scheme to different applications 

Stability analysis and checking performed to the scheme in the general cKdV equations give the ability 
of using this scheme to solve other problems for which analytical solutions have not been found. We first 
consider the multi-soliton solution decay of the localized initial condition for the single KdV equation of 
HS system (figure 4. a) then for the complete HS system (figure 4.b). In both we use the initial condition 
from formula (14) but with 10 times the width and twice the amplitude. 





First KdV 
i 

; i 

2 5 i 


esq li git 

"7 - 
6 


ion 


of HS system 

t=o 

t=Ssec 






:: t ! i 

1 1 1 1 u 


5 - 












■s AX 2 








-60 


-40 -20 


o — 

c 


5 X 


20 40 


60 



Fig. 4. A The multi-soliton decaying of the isolated (first) KdV of HS system (13. a) 




Fig. 4.b The multi-soliton decaying of the complete HS system (13). 

The second mode affects (interacts) the first one which results in the right direction soliton-like "tail" as 
shown in the figure 4.b. We estimate, using the conservation low (derived below), the expected number of 
solitons that already obtained in the numerical solution as 

0i x (13.1) - 20 2 x (13.2) =>- 
I [-50? - 01} + £[-.59f - .250^ + .1250^ - 6 2 6 2xx + .50 2 2 J = 
6i and 6 2 — > as x — > ±oo 
ioo°° (-5^1 — 6l)dx = const. 

Next we go to the solution of nonintegrable HS system. The integrable HS system (13.a,b) may be shifted 
to "slightly" nonintegrable one by small change of the dispersion constant of the first equation to have the 
new nonintegrable HS system 

(0O t - 0.2 (0O 3 , - 1-5 (00* (0i) + 3 (6 2 ) x (0 2 ) = (15.a) 



(0 2 ) t + 0.5 (6 2 ) 3x + 1.5 (6 2 ) x (9 1 ) = 0. 



(15.b) 



Using our scheme with initial condition from (14) we find that the scheme works satisfactory (in the sence 
of convergence) even for nonintegrable HS system as shown from Figure 5 below. The solution looks like 
a soliton one for small time. 



Firstmode, t=0 

t=6sec(integrable) © 

— — t=6see( non-integ rable^ 



Second mode, t=0 

t=6sec(integrable) 6 - 

— t=6sec(non-integrable)5 - 




Fig. (5) The numerical solution of integrable and slightly nonintegrable HS system. 
Also the solution using a non-smooth initial condition for HS system (13) is shown in figure 6 below 
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Fig.6 The numerical solution of Non-smooth initial condition for HS system (13). 
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